eval '(exit $?0)' && eval 'exec perl -S $0 ${1+"$@"}' && eval 'exec perl -S $0 $argv:q' if 0;
#;-*- Perl -*-

# 13-05-2005

# Makes a linear interpolation between to POSCARs with N points

use FindBin qw($Bin) ;
use lib "$Bin" ;
use Vasp ;

(@ARGV<3) || die "USAGE: getcif.pl <POSCAR POTCAR>\n" ;
# Defaults
$poscarfilename="CONTCAR";
$potcar="POTCAR";
# Read input variables
if(@ARGV>0){
  $poscarfilename=$ARGV[0];
}
if(@ARGV>1){
  $poscarfilename=$ARGV[0];
  $potcar=$ARGV[1];
}
(-e "$poscarfilename") || die (" $poscarfilename does not exist\n");
(-e "$potcar") || die (" $potcar does not exist\n");
($coordinates,$basis,$lattice,$num_atoms,$total_atoms,$selectiveflag,$selective,$description)
 = read_poscar($poscarfilename);
$coms=`grep 'TITEL' $potcar`;
@coms=split /\n/, $coms;
$num_coms=@coms;
(@{$num_atoms}==$num_coms) || die (" POSCAR and POTCAR files should have the same types of atoms\n");
for($i=0;$i<$total_atoms;$i++){
  for($j=0;$j<3;$j++){
  if($coordinates->[$i][$j]<0) {
    $coordinates->[$i][$j]+=1;}
 if($coordinates->[$i][$j]>=1) {
    $coordinates->[$i][$j]-=1;  }}}
open OUT, ">out.cif";
printf"data_\n_audit_creation_method    ''\n";
printf OUT "data_\n_audit_creation_method    ''\n";
$lengtha=sqrt($basis->[0][0]*$basis->[0][0]+$basis->[1][0]*$basis->[1][0]+$basis->[2][0]*$basis->[2][0]);
$lengthb=sqrt($basis->[0][1]*$basis->[0][1]+$basis->[1][1]*$basis->[1][1]+$basis->[2][1]*$basis->[2][1]);
$lengthc=sqrt($basis->[0][2]*$basis->[0][2]+$basis->[1][2]*$basis->[1][2]+$basis->[2][2]*$basis->[2][2]);
$cos_alpha=($basis->[0][1]*$basis->[0][2]+$basis->[1][1]*$basis->[1][2]+$basis->[2][1]*$basis->[2][2])/$lengthb/$lengthc;
$sin_alpha=sqrt(1-$cos_alpha*$cos_alpha);
$alpha=atan2($sin_alpha,$cos_alpha)*180/3.141592653589793;
$cos_beta=($basis->[0][0]*$basis->[0][2]+$basis->[1][0]*$basis->[1][2]+$basis->[2][0]*$basis->[2][2])/$lengtha/$lengthc;
$sin_beta=sqrt(1-$cos_beta*$cos_beta);
$beta=atan2($sin_beta,$cos_beta)*180/3.141592653589793;
$cos_gamma=($basis->[0][0]*$basis->[0][1]+$basis->[1][0]*$basis->[1][1]+$basis->[2][0]*$basis->[2][1])/$lengtha/$lengthb;
$sin_gamma=sqrt(1-$cos_gamma*$cos_gamma);
$gamma=atan2($sin_gamma,$cos_gamma)*180/3.141592653589793;
printf"_cell_length_a%20.6f\n",$lengtha;
printf OUT "_cell_length_a%20.6f\n",$lengtha;
printf"_cell_length_b%20.6f\n",$lengthb;
printf OUT "_cell_length_b%20.6f\n",$lengthb;
printf"_cell_length_c%20.6f\n",$lengthc;
printf OUT "_cell_length_c%20.6f\n",$lengthc;
printf"_cell_angle_alpha     %.4f\n",$alpha;
printf OUT "_cell_angle_alpha     $.4f\n",$alpha;
printf"_cell_angle_beta      %.4f\n",$beta;
printf OUT "_cell_angle_beta      $.4f\n",$beta;
printf"_cell_angle_gamma     %.4f\n",$gamma;
printf OUT "_cell_angle_gamma     %.4f\n",$gamma;
printf"_symmetry_space_group_name_H-M  'T'\n";
printf OUT "_symmetry_space_group_name_H-M  'T'\n";
printf"loop_\n";
printf OUT "loop_\n";
printf"_atom_site_type_symbol\n_atom_site_label\n";
printf OUT "_atom_site_type_symbol\n_atom_site_label\n";
printf"_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n";
printf OUT "_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n";
$index=0;
for($j=0;$j<@{$num_atoms};$j++){  
     $line=$coms[$j];
     chomp($line) ;
     $line=~s/^\s+//g;
     @line=split /\s+/,$line;
     $c[$j]=$line[3];
    if($j != 0){$num_atoms->[$j]+=$num_atoms->[$j-1];}  
for($m=$index;$m<$num_atoms->[$j];$m++){
printf "%5s%5s%d\t",$c[$j],$c[$j],$m+1;
printf OUT  "%5s%5s%d\t",$c[$j],$c[$j],$m+1;
printf "%12.8f%12.8f%12.8f\n",$coordinates->[$m][0],$coordinates->[$m][1],$coordinates->[$m][2];
printf OUT  "%12.8f%12.8f%12.8f\n",$coordinates->[$m][0],$coordinates->[$m][1],$coordinates->[$m][2];
}
$index=$num_atoms->[$j];
}


